Novel self-assembled morphologies from isotropic interactions 
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We present results from particle simulations with isotropic medium range interactions in two 
dimensions. At low temperature novel types of aggregated structures appear. We show that these 
structures can be explained by spontaneous symmetry breaking in analytic solutions to an adaptation 
of the spherical spin model. We predict the critical particle number where the symmetry breaking 
occurs and show that the resulting phase diagram agrees well with results from particle simulations. 
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Understanding the principles behind spontaneous for- 
mation of structured morphologies is of interest both as 
a fundamental scientific question and in engineering ap- 
plications where the possibility of using self-assembly to 
produce novel materials provides a compelling comple- 
ment to traditional blue-printed fabrication [l|, |2[ . Con- 
sequently there is growing interest in exploring interac- 
tions that can facilitate self-fabrication of materials with 
novel properties |3|— lSl] • The more general question of what 
structures are possible to self-assemble from a given class 
of interactions has however not been addressed, with 
some notable exceptions e.g. simulation based analysis 
of polyhedral packing [6j. 

To examine the possibilities of self-assembly from a 
theoretical angle, in this Letter we consider systems with 
pair- wise isotropic interactions, frequently used as coarse 
grained models of more complex meso-scale systems such 
as colloidal systems 0, particles in an ambient fluid (HQ? 
and spin glasses [10j. We show that typical aggregates 
appearing as low temperature configurations in parti- 
cle simulations with randomly generated medium range 
isotropic interactions can be predicted analytically us- 
ing an adaptation of the spherical spin model [11]. The 
morphologies, many of them novel and surprisingly com- 
plex, can be systematically classified by their sponta- 
neous breaking of the rotational symmetry. 

To formulate a solvable model of a particle system we 
start by considering a lattice spin system with Hamilto- 
nian on the form 



H 



(i) 



where si G {0,1}, and Si = 1 represents a particle at 
lattice site i and = represents vacuum. The total 
number of particles is set by the normalization ^ Si = 
TV. The interaction U^ is an effective isotropic potential 
with a hard core repulsion corresponding to the lattice 
spacing. In a spin glass metal the interaction is mediated 
by a polarization of the Fermi sea [l2[ while in a colloidal 
system it could involve a surface polymer induced steric 
hindrance competing with a depletion attraction [l3| . 
The discrete model described by Eq. (pQ) can equiva- 



lent ly be formulated as a continuous one (si G R) with 
auxiliary constraints, Yli s 7 ^ = N Vra. To make the 
model analytically tractable we relax the constraints to 
include only the first two moments. The result is the 
spherical spin model including an external field, with 
Hamiltonian H s = J^. UijSiSj + h^2 i Si. A neces- 
sary condition for an energy minimum in this model is 
J2j(Uij — ^Sij)sj = h/2, where 5 is the identity matrix, 
and A comes from the constraint on the second moment. 
Translational invariance implies that S{ = c is a solution. 
Non-constant solutions Si = v\ + c exist if A is an eigen- 
value of the matrix U corresponding to the eigenvector 
v. From the constraints it follows that c = p, where p 



is the density, and V 1 ^ 



(1 — p)p, where V is 



the number of lattice sites or the volume. The ground 
states are defined by the eigenvector (s) corresponding to 
the lowest eigenvalue of U. 

It can be shown, e.g. using translational invariance, 
that all isotropic interaction matrices have the Fourier 
normal modes as eigenvectors. This also follows from 
the observation that any isotropic interaction matrix Uij 
can be expressed as a linear combination of matrices on 
the form A^A™, /, m = 0,1,..., where the two compo- 
nents are the discrete Laplace operator defined as usual 
A + or along the diagonals A x (note that A + and A x 
commute). The Fourier modes are eigenfunctions of the 
Laplace operator so the two arguments lead to the same 
conclusion. However, the latter argument points to a 
subtlety: the eigenfunctions of the Laplacian depend on 
the boundary conditions, where the Fourier harmonics 
result from periodic boundaries. Requiring that the func- 
tions converge to zero at infinity instead results in cylin- 
drical (spherical etc.) harmonic eigenfunctions with an- 
gular modulations localized around a nucleation point. 
These localized configurations bears similarity to topo- 
logical defects, such as the Belavin-Polyakov monopole, 
that are important for the low temperature behavior of 
e.g. Heisenberg ferromagnets [14]. 

Our central result is that, when there are too few par- 
ticles to occupy a translational invariant ground state, 
the eigenfunctions of Bessel type J u; (27r/cr) cos (ujO) de- 
termine the low temperature behavior of self-assembling 
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FIG. 1. Predicted morphological alphabet for aggregating potentials generated by equation (3]) (with uj < 8). The morphologies 
shown are energetically preferable for a large fraction of random potentials. Red background indicates a strong signal while 
blue background signifies weaker signal, a denotes the number of active radial oscillations of the ground frequency after the 
threshold mapping. 



particle systems in two dimensions. The eigenvalues, i.e. 
the energies, are given by the radial Fourier transform of 
the potential as 15] 16| 



E(k) = 2tt [ rdr U(r)J {2nkr), 
Jo 



(2) 



so that the wavenumber k = arg min fe E(k) defines an 
eigenfunction with minimal energy. The energy spec- 
trum is degenerate since E(k) is independent of lj, reflect- 
ing the rotational invariance around the nucleation point. 
However, this degeneracy is broken by the non-linearity 
in the mapping from the spherical model to the discrete 
lattice system, where a single base frequency lj together 
with its overtones (similar to the harmonic overtones of a 
square wave) are energetically favored. The spontaneous 
symmetry breaking of the ground state, from 0(2) to one 
of its isotropy subgroups , is equivalent to the behav- 
ior in a bifurcation problem 17|, LL8| . The ground states 



of the spherical model of relevance to the self-assembly 
problem in (pp) are therefore on the form 



a n J noJ (27rtvr) cos(nuj0) + cJo(27r/c m ) 



(3) 



n=0 



in the limit k m — >• 0, which turns the last term (which 
we refer to as the mass-builder) into the translational in- 
variant constant c in the minimization of H s . However, 
the k m —> limit is only relevant for infinite structures 
(with non-zero global density). For an aggregated struc- 
ture with finite mass a small but non-zero k m is needed to 
localize the solution (resulting in a non-zero local density 
but a global density of zero). 

As we will see later, the sum in Eq. ((3]) is dominated 
by the first few terms and in practice most structures 
are adequately described by the zeroth and first term. 



Notable exceptions are finite lattices, see e.g. Fig. QJ, 
which grow by successively including higher terms. 

To connect the above results to the particle model (pQ), 
we map the continuous ground states (si G M) to particle 
configurations (si G {0, 1}) by applying a step function 
threshold. The energy of the resulting discrete configu- 
rations depend on how the mapping distorts the power 
spectrum relative to the energy spectrum ([2]). While the 
main weight of the discrete configurations typically re- 
mains localized to k, contributions from overtones that 
appears in the discrete system can have a large effect 
on the energy. These effects are difficult to quantify an- 
alytically. Instead we generate possible candidates for 
ground states using Eq. (|3|), i.e. linear combinations of 
Bessel functions with angular frequencies 0, cj, 2uj and 
3uj with k = 1, as well as a mass-building Bessel func- 
tion with uj = and k m chosen such that the first zero 
of the mass-builder is slightly larger than the ath zero of 
Juj for a = 1,2,3. The configurations were mapped to 
binary valued configurations on a lattice with a thresh- 
old that defines different masses. The energy of the re- 
sulting configurations were calculated for 1000 random 
first-order spline potentials. We used aggregating poten- 
tials, i.e. potentials with a global minimum in the energy 
spectrum at k = 0. In addition their arguments where 
scaled so that the interior minimum in the energy spec- 
trum resided at k = 1 as demonstrated in Fig. [2^ and 
were attractive at the lattice distance to ensure correct 
resolution of the lattice. For each potential and mass we 
recorded the configuration with the lowest energy, which 
resulted in a limited number of favored morphologies, 
shown in Fig. [TJ The method predicts many novel struc- 
tures (e.g. Fig. QJi,c,d,g,h,m), as well as simpler ones like 
disks and localized lattices (e.g. FigJTj3,e,l,q). The latter 
are observed on the atomic scale [l9j], while the former 
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FIG. 2. Showing self- assembled morphologies for a single potential, (a) A typical energy spectrum of a random potential. The 
wavelength of the morphologies is determined by the interior minimum as of the spectrum when the number of particles N is 
less than iV c , the number required to build a disk with radius corresponding to ko, see Eq. J4|). (b)-(g) Transitions between 
different symmetry groups for the potential of (a) as the particle number is decreased. (Top) Configurations from annealed 
particle simulations, (middle) corresponding states of the spherical model with contours added at a level giving correct number 
of particles, and (bottom) Bessel spectra (arbitrary units) of the particle configurations expanded on the form of ©, summed 
over uo. The first maximum corresponds to the mass-builder /c m , the second should be compared to as = 1 predicted from the 
spectrum in (a), and the rest are overtones. 
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FIG. 3. Examples of configurations from particle simula- 
tions with randomly generated potentials. All configurations 
except (b) are well represented by functions on the form Q 
with the first four terms included, most requiring only one or 
two. (a) Morphologies predicted as common by our method, 
see Fig. [T] (b) An example of a hierarchical structure, not di- 
rectly describable in our theory, where each part of the usual 
Bessel structure serves as the nucleation point for a separate 
Bessel structure, (c) Morphologies that, while being simply 
representable by Eq. ([3]), are not predicted in Fig. U 



morphologies are only expected to appear at higher par- 
ticle numbers. 

To summarize the predictions of the adapted spherical 
spin model, we expect a particle system with an aggregat- 
ing potential to show (I) morphologies from a limited al- 
phabet (Fig. [I]) with (II) a single dominant wavelength n 
determined by the minimum of the spectrum and (III) an 
uo degeneracy where a single potential can self-assemble 
into many different morphologies, depending on external 
parameters. 

To test these predictions we performed off-lattice simu- 
lations of particle systems. We constructed 1200 random 
interaction potentials (piecewise constant and 3rd order 



splines with the same restrictions on the spectrum as 
above) and did Monte Carlo annealing at different parti- 
cle numbers. Examples of the resulting particle configu- 
rations are shown in Fig. [2]|3] and are in good agreement 
with (I). Of the simulated particle systems approximately 
86% (Fig. [3^) annealed to predicted morphologies shown 
in Fig. Q] and 9% (Fig. [3fc) annealed to morphologies de- 
scribed by Eq. ([3]) but absent in Fig. [T] due to the limited 
parameter range considered (e.g. uo < 8). The latter were 
in general similar to those in Fig. [TJ The frequencies of 
which predicted configurations where observed was also 
largely consistent with the strength of the signal of the 
corresponding morphology in Fig. [TJ The remaining 5% 
of the morphologies exhibited hierarchies of nested Bessel 
like structures, for an example see Fig. [3b. 

By expressing the observed particle configurations in 
the modulated Bessel basis we confirm that (II) the dom- 
inant wavenumber is accurately predicted by the mini- 
mum in the energy spectrum of the potential. Examples 
of this are shown in Fig. EJ>g, where we note that while 
the point symmetry of the preferred configuration for a 
single potential changes with varying particle number, 
the dominant (non mass-building) active wavelength is 
approximately constant and equal to the predicted value 
n. This switching of symmetry group is observed fre- 
quently in the simulations and shows that (III) the uo 
degeneracy is not just a mathematical curiosity in our 
model. For a given potential the particle number sets 
the morphology, an important fact that could facilitate 
the use of standard techniques such as density gradient 
centrifugation [2o| for differential selection of morpholo- 
gies 

The spectral analysis also allow us to analytically pre- 
dict a phase diagram for aggregating potentials. The 
spectrum shown in Fig. [2^ is typical for the potentials 
we consider in that it is oscillatory and has a global min- 
imum at k = 0. In a spin system a global minimum at 
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FIG. 4. Decreasing the number of particles causes a tran- 
sition from disc-like structures at the macro-scale to more 
complex morphologies at the meso-scale. The figure shows 
the phase diagram with the predicted (full line) and mea- 
sured (markers) critical number of particles N c , adjusted for 
packing fraction /, versus ko calculated (see Fig. [2Jl) from 
respective spectra for 33 random potentials. (Inset) Four ex- 
amples of potentials used. 

k = implies a ferromagnetic ground state. For a par- 
ticle system it implies that in the high particle limit we 
can use large scale interface minimization to argue that 
the ground state is dominated by the mass-builder, typi- 
cally a closely packed disc with a surface that is smooth 
or has wavelike indentations (similar to Fig. [2t). For a 
given number of particles there is a maximal size of such 
a disc and for low particle numbers this effectively ex- 
cludes the small k part of the energy spectrum in the 
minimization that determines k. This causes a transi- 
tion as the particle number is decreased, from a disc to 
the more complex Bessel based morphologies we observe. 
This transition happens when the energy of the disc raises 
above the energy of the lowest interior minimum E(k). 
We denote this point fco, as illustrated in Fig. [2^. The 
critical particle number N c can be predicted through 



where / is the packing fraction, a is the radius of the 
particles, and b « 0.293 is a numerical constant corre- 
sponding to the optimal k to describe a disc of radius 1, 
found through arg max fe J* Jq (2irkr) rdr / 1 1 Jo (fc) 1 1 \ . The 
phase diagram predicted by Eq. (j4j) is shown in Fig. |4] 
together with simulation results. 

We conclude that isotropic pairwise additive potentials 
can give rise to complex morphologies appearing between 
the atomic and macroscopic scale and that the frequently 
occurring structures forms a limited morphological al- 
phabet. The patterns are generically constructed as dis- 
cretized linear combinations of a few Bessel functions, 



which can be understood and analytically predicted from 
an adaptation of the spherical spin model. Further, we 
analytically calculate the phase diagram showing where 
these patterns occur for different potentials. The accu- 
racy of the predictions is surprising considering the com- 
plexity involved in determining the ground states in a 
many particle system. The methodology we present ap- 
plies to the entire class of isotropic interaction poten- 
tials and provides new theoretical understanding of self- 
assembly processes. 
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